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Abstract. We solve the integral equation describing the propagation of light in an isothermal plane-parallel atmosphere of 
optical thickness t*, adopting a uniform thermalization parameter e. Fhe solution given by the ALI method, widely used in the 
field of stellar atmospheres modelling, is compared to the exact solution. Graphs are given that illustrate the accuracy of the 
ALI solution as a function of the parameters e, t* and optical depth variable r. 
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1. Introduction 

The solution of the radiative transfer equation (RTE) is at the 
heart of the stellar atmospheres modelling, since this equation 
has to be solved typically thousands of times in order to con- 
struct a realistic model. It is thus crucial to get a clear idea of the 
accuracy with which the RTE is solved, and the effect it has on 
the determination of the main physical quantities of the model: 
populations, electron density, temperature, etc. In this article, 
we focus on the first point checking the ALI method for solving 
the integral form of the RTE, since this method is nowadays at 
the basis of most numerical schemes used to determine the radi- 
ation field in stellar atmospheres. We recall that "ALI" means 
Accelerated (or Approximate) Lambda Iteration, the Lambda 
operator being defined by Eqs. below for the scattering 

law we adopt here. The ALI code used in this paper is a com- 
bination of an accelerated iterative method (with a diagonal A- 
operator) and a formal solver based on parabolic short charac- 
teristics. Recent reviews on this approach are Paletou (2001), 
Hubeny (2003) and section 3 of Trujillo Bueno (2003). 

The accuracy of our ALI code is tested while applied to a 
well-known problem consisting of a homogeneous, isothermal 
slab with isotropic and monochromatic light scattering (Sec.|2). 
Indeed, this idealized problem can be solved exactly, which al- 
lows for a direct comparison with the solution given by the 
ALI method. This problem is very simple on physical grounds 



but implies analytical and numerical calculations that are far 
from trivial. It contains the seeds of most of the difficulties met 
when solving the RTE in a thick, highly scattering medium. 
It thus provides an excellent test for numerical codes since 
very accurate analytical solutions are available (Sec. [3jl. After 
a brief description of our ALI code, we move to the numerical 
tests in Sec. @J which is the main part of this paper. The link 
with previous studies on the subject (Trujillo Bueno & Fabiani 
Bendicho 1995, Trujillo Bueno & Manso Sainz 1999) is finally 
commented in Sec. [5] 



2. The standard radiative transfer problem 

This problem consists in solving the RTE in a homogeneous 
plane-parallel atmosphere of optical thickness r* > (pos- 
sibly infinite); light scattering is assumed to be isotropic and 
monochromatic. It is furthermore supposed that the matter is 
in local thermodynamical equilibrium with uniform tempera- 
ture T through the atmosphere. The thermal source function at 
any frequency is then eB(T), where e is the (spatially invari- 
ant) photon destruction probability per scattering and B(T) the 
Planck function at temperature T (frequency dependence is not 
mentioned). In the absence of any external source of radiation, 
this problem reduces to solving the following integral equation 
for the source function S (Mihalas 1978): 
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S(T) = efl(r) + (l-6)(AS)(T), 



(1) 
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where the A-operator for isotropic and monochromatic scatter- 
ing is 

(A5)(r) = iJ E 1 (\T-T'\)S(j / )dT'. (2) 

Here, E\ is the first exponential integral function as defined by 
d/i 



£i(t) 



exp(-r//i) — 



(t > 0) . 



(3) 



We remind the reader that Eq. Q models the multiple scatter- 
ing of photons of frequency v assuming that 1) the scattering 
is monochromatic (or coherent) if v belongs to a continuum, 2) 
the line profile is rectangular (Milne profile) if v belongs to a 
spectral line (see, e.g., Ivanov 1973, p. 57). 

The solution to problem Q is S(r) = S(e,T*,T)B(T), 
where S (e, r*, r) satisfies the integral equation 



S(e,r*,T) = e + (l -e)(AS)(e,T*,-r) 



(4) 



depending on parameters e and r*. Note that this function is 
symmetrical about the r-mid-plane: S(e, t*,t) = S(e, t*,t* - 
r). 

This equation is the integral formulation of the RTE in our 
model; it specifies the standard radiative transfer problem we 
intend to solve analytically (Sec.[3jl and numerically (Sec.|4j. 

3. Analytical solution of the standard problem 

There are many analytical methods for solving the integral 
equation The classical approach, recently reviewed by 
Chevallier & Rutily (2003, hereafter Paper I), involves the basic 
auxiliary functions of radiative transfer theory in plane-parallel 
geometry, namely the //-function for a semi-infinite space, and 
the X- and T-functions for a finite slab (Chandrasekhar 1960). 
The //-function depends on the parameter e and on an angular 
variable taken as positive hereafter. In addition the X- and 
F-functions depend on t*, and we have X{e,T*,fj) — > H(e,/S) 
and Y{e,T*,fi) — > as r* — > +oo. 

The zero-order moments of the functions H, X, and Y yield 
the surface values of the solution S to @. The moment of the 
//-function is defined and given by 



ao(e) 



Jo 



1 + yfe 



(5) 



and those of the X- and T-functions defined as 

a (e,T*)=f X(e,T*, M )dn, A,(e,r*)=f T(e,T»d^ (6) 
Jo Jo 



are related by 



1 - — r-aoO.T*) 



1-e 



: A)(e,r*) 



(7) 



There is no exact expression of these moments. 

In a semi-infinite atmosphere, the surface value of the solu- 
tion S to ® is 



S(e,0) = 1 



1 



-a (e) = ^ 



(8) 
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Fig. 1. Source function S (e, t*, t) in a slab with e = 10 
2000. Right and top axes display the linear plot with infinite 
first derivative at r = 0. Black dots show the thermalization 
depth l/k(e) on both scales, the definition of k(e) being given 
in Paper I. 



and it is 



S(e,T*,0) = 1 



-[a (e,T*)+A)(e,r*)] 



(9) 



in a finite slab. As S is symmetrical about the r-mid-plane, 
S(e, t*,t*) = S(e, t*,0). These relations were first derived by 
Sobolev (1957, 1958). The former result is the famous 
law" for semi-infinite media. The latter one is less known; it 
requires a table of moments (ao, /3o) for numerical applications. 
Such tables are available in the literature: see references in Van 
de Hulst (1980, p. 225-227). Very accurate surface values of 
the S -function can also be found in Paper I. 

The calculation of the function S within the slab is dis- 
cussed in detail in Paper I, which contains ten-figure tables 
of S(e,T*,T) for (e,r*) = (0.5,2), (10" 2 ,20), (10" 4 ,2000) and 
(10~ 8 ,2 x 10 s ). In a half-space, the internal solution is known 
since the end of the 50's and it can be expressed in closed-form 
in terms of the //-function. In a finite slab, the solution involves 
two non-classical auxiliary functions £ + and that are implic- 
itly defined by Fredholm integral equations over [0,1]. These 
equations can be solved very accurately, so that the solution in 
a finite slab is nearly as accurate as in a half-space. The accu- 
racy is estimated at better than 10~ 10 for any value of e, t* and 
t, which means that the solution given in Paper I can safely be 
used as an accuracy test of the ALI code. 

The general behavior of the S -function is shown in Fig.^ 
which illustrates the Table 3 of Paper I (e = 10~ 4 and r* = 
2000). It can be seen that the solution S tends to 1 for large val- 
ues of t and that it drops when r is close to the thermalization 
depth 1 /k(e) « 58 for e = 10 4 , where k(e) is defined in Paper 
I. It tends steeply to the surface value S (0) as r tends to 0, with 
an infinite derivative at 0. It is regrettable that the generally 
adopted logarithmic scale in r obscures this essential last point, 
as seen when comparing the solid and dashed curves of Fig.^ 
The explanation lies in the fact that dS / <9(log t) = t dS jdr — > 
even if dS/dr ~ E\(j) — > +oo. 
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4. Comparison with ALI numerical solutions 

In this section we compare in detail the analytical solution 
described in the previous section to the one given by our 
ALI code. This code uses a diagonal approximate A-operator 
(Olson et al. 1986). At each iteration, a formal solver has to 
be used in order to calculate the transform of the source func- 
tion by the A-operator. Inserting the definition of the E\- 
function into Eq. (0 and inverting the order of integrations, the 
so-called formal solution to the RTE is first calculated 



1 



— S(r')e X p[(T-T')/yu] dr' 

y" Jo 
S(t) 



if - 1 < n < 0, 
if H = 0, 

- f S(t') exp [-(t' -t) In] At' if < fi < +1. 

(10) 



Then the A-transform of the source function is derived, since it 
is here the associated mean intensity 



(AS)(t) 



2 J-i 



I(j,jd)diJ.. 



(ID 



The formal solution (II Ui is calculated following the method 
of short characteristics whose basic elements can be found in 
Olson & Kunasz (1987) and Kunasz & Auer (1988). It was 
further improved by the implementation of monotonic inter- 
polation for multi-dimensional applications (Auer & Paletou 
1994) and by Fabiani Bendicho & Trujillo Bueno (1999) for 
three-dimensional applications with horizontal periodic bound- 
ary conditions. In the present paper, we used parabolic short 
characteristics. The //-integration in (1111 is performed with the 
help of a Gaussian quadrature. 

A numerical acceleration scheme is used so as to im- 
prove the rate of convergence of ALI: this is the so-called Ng- 
acceleration introduced in the field of radiative transfer by Auer 
(1987, 1991; see also Rybicki & Hummer 1991). 

We have calculated the relative error 



d(e, t*, t) = 



SAU(€, T*,T)-S(e,T*,T) 



S(e, t*,t) 



(12) 



at various optical depths, where S (e, r*, t) is the analytical so- 
lution of Sec. [3]and S ALi(e, t*, t) is the solution given by the 
ALI code. This error corresponds to the "true error" defined by 
Auer, Fabiani Bendicho & Trujillo Bueno (1994), who used a 
finer grid to calculate S (e, t*, t). 

We introduce also the maximum value of d{e,T*,T) when 
the r-variable covers the domain [0, r*], viz. 



dyi{e, T*) = max d(e,T*,r). 

0<T<T* 



(13) 



Of course d and dyi depend on the number of iterations 
N performed by the ALI code during each run. Finally we 
define N c as the number of iterations used to reach conver- 
gence, which is the smallest value of N satisfying the condition 
|1 - d M (N)/d M (+oo)\ < s c , where d M (+oc) = d M (N = 10000) 
and e c is arbitrarily set to 0.01 in the present paper. 

The slab optical depth is discretized using a logarithmic 
grid, symmetric with respect to the mid-plane, with « T points 
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Fig. 2. Relative error a!(e, t*,t) of the ALI code for n T - 9, 
«// = 5, r m = 10 4 and N = 1000. Black dots show the ther- 
malization depth 1 /k(e) for each case. 



per decade, including the r = point, the next point denoted 
by T m , and the last point t - t*/2. The angular integration 
in Eq. il 1> is performed with a symmetric grid containing 
Gauss-Legendre points in [0,1]. There is no frequency integra- 
tion since light scattering has been supposed monochromatic. 
In most of our calculations, we chose the values r m = 10~ 4 , 
n T = 9, and — 5 (i.e., values quite often adopted for stellar 
atmospheres modelling). Some values of (e, r*) may be (0.01, 
20) for a continuum, (10~ 4 , 2000) for an "average" spectral line 
and (10~ 8 , 2 x 10 8 ) for a strong spectral line. 

The quantities of interest are the maximum relative error du 
and the number of iterations to reach convergence A^ c , which 
depend on e, t* and numerical parameters r m , n T , and N. 
We first study the variation of d(e, t*, t) with t. Then, we study 
the influence of e, r*, n T on du and A^ c , for given r m , and for 
N = N c . 

4.1. The influence ofr and N 

Figure|5]shows the variation of the relative error t — » die, t* , t) 
for the three selected values of e and r*. It can be seen that 
the accuracy (i.e. maximum relative error) of our ALI code is 
about 5 x 10~ 3 for the three cases studied here. Relative error is 
close to this accuracy when t is smaller than the thermalization 
depth l/k(e) of the atmosphere (black dots on the curves), and 
significantly improves beyond (up to 10~ 8 ). In photon mean 
free path units, the thermalization depth is 6, 58 and 5774 for 
e = 10~ 2 , 10~ 4 and 10~ 8 respectively. Note that the surface rel- 
ative error is a good estimator of the accuracy in spectral lines, 
but not in the continuum. 

The iterative algorithm was stopped after = 1000 itera- 
tions. Figure |3] shows that this number ensures convergence of 
the ALI code, the convergence being slower when e — > and 
t* — > +oo. Irregular steps in these curves are due to the Ng 
acceleration process, here operated every four iterations. 

4.2. The influence of r m 

We point out that du is improved when r m goes to 0, up to a 
given value where the accuracy is constant. Including the t = 
point in the grid and choosing r m < 10~ 2 , the best accuracy 
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Fig. 3. Evolution of the maximum relative error dyi with itera- 
tion number N for n T - 9, «„ = 5 and r m = 1CT 4 . 



is warranted. Excluding the t = point from the grid has no 
influence on accuracy if we choose r m < 10~ 4 . The standard 
choice r m = 10~ 4 is thus correct, and this value will be adopted 
hereafter. 



4.3. The influence ofn T andn M 

In Fig. 0] are shown the variations of du in an average line as 
a function of n T and n^. The maximum relative error a?M de- 
creases with increasing number n T of T-grid points, and it is 
sensitive to the choice of the number n M of angular grid points 
up to an optimal value «J,° pt) ; the latter is defined as the smallest 
value of rifj for which the condition 1 1 - du(/i/^/d\i(64)\ < 0.01 
holds. The accuracy does not increase with a finer yU-grid. We 
note that n|, opt) < n T and that we have a linear dependence of 
this optimal value on n T : «^ opt) = 0.8n T + 2.7 (dashed curve). 
This fit is still valid for strong lines. 

As seen in Fig. [5] (same as Fig.|4]for the continuum), the 
fit for lines cannot be applied to the continuum, for which 
nf pl) > n T . It is still possible to define and calculate an optimal 
value for n T < 18, using the relation n ( ° pi) = 1.7n T + 1.3 (dashed 
curve). It appears that our ALI code is more demanding in an- 
gular resolution when solving the problem @ in a continuum 
than in a line. 

The results of Figs. 0] and [5] are detailed in Fig. [6] for the 
three chosen values of (e, f*) and = 64. We remark that the 
accuracy improves with n T for each couple (e, t*), more signif- 
icantly in the continuum than in lines. Figure0gives the num- 
ber of iterations N c used to reach convergence for e c = 10~ 2 and 
rt/j = 64. The number N c appreciably increases with n T in the 
lines: it is indeed well known that the rate of convergence of the 
one-point ALI iterative scheme drops for an increasing refine- 
ment of the spatial grid (Olson et al. 1986); however improve- 
ments were already proposed (e.g., Trujillo Bueno & Fabiani 
Bendicho 1995) in order to increase significantly the rate of 
convergence of ALI-based methods. 

4.4. The influence of e and t* 

The maximum relative error du and number of iterations A^ c 
are shown in Figs.|S]and|5]for an extended range of (e, r*) after 




Fig. 4. Maximum relative error du in an average line (t* = 
2000, e = 10~ 4 ) as a function of n T and n^. The dashed curve 
represents a linear fit n|f pt) = 0.8« T + 2.7. 




Fig. 5. Maximum relative error a?M in the continuum (r* = 20, 
e = 0.01) as a function of n T and n M . Dashed curve represents a 
linear fit n° pt> = 1 .ln T + 1 .3 for n T < 18. 



the ALI code has converged (N - 10000, n T — 9 and — 5 
are fixed here). 

As seen in Fig. [8] the accuracy hardly changes as e — > 
and t* — > +oo, but the number of iterations needed to achieve 
convergence increases substantially (see Fig.|9j- When e > 0.1 
the accuracy no longer depends on values of r*. The compari- 
son of Figs.[6]and[8]leads to a disagreement since the parameter 
rifj is set to different values, 64 and 5 respectively. 

In Fig. [9] we have plotted the parameter A^ c as a function of 
e and r*. When e — > and t* — > +00, we note a slowing down 
of the convergence (already seen in Fig.|3}- 

5. Comments on previous studies 

Now we compare our results for the monochromatic scatter- 
ing problem with those published by Trujillo Bueno & Fabiani 
Bendicho (1995) and Trujillo Bueno & Manso Sainz (1999). 
Although these two papers concern mainly the development of 
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new iterative methods for radiative transfer applications (for the 
unpolarized and polarized cases respectively) they give some 
information on the accuracy of the numerical solutions ob- 
tained for spatial grids of increasing resolution. 

In Table [lj good agreement is found between our values of 
N c , d(e, t*, 0) and those given by these authors; our surface rel- 
ative error die, f, 0) corresponds to their surface true error T e . 
The observed small discrepancies are possibly due to the dif- 
ferent scattering laws adopted, leading to different A-operators. 




Fig. 9. Number of iterations iV c (e c = 10 2 ) as a function of e 
and r* (n T = 9, = 5). 
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Fig. 10. Relative difference Si = 1 1 - S (e, t*, 0)/ s/e\ as a func- 
tion of k(e)T* . Dots represent values of e e [10~ 12 , 1] and 
t* G [0.2, 2 X 10 8 ], and the solid curve is exp (-£(e)r*). 

In Fig.^|it is shown how far the semi-infinite exact result 
S(e, 0) = Vir agrees with the finite one. This comparison is 
useful since many authors use the -\fe-law as a check for their 
calculations in thick slabs. We plot the relative difference 5i = 
| 1 - 5(6,^,0)/ yfe\ as a function of k(e)T*, and the quantity 
exp (—k(ejr*) which characterizes the validity of the -\fe-law 
(solid curve). For k(e)T* > 100, the accuracy limit of our code 
is reached, which explains that the solid curve no longer fits the 
dots. This law is very well satisfied in lines (fe(e)r* * 34 in an 
average line) but not enough in the continuum (fc(e)r* m 3.3). 
We conclude that the -\Ze-law can be used as a test for the ALI 
code when k(e)r* > 10, since then is an approximation to 
the surface value with an accuracy better than 10~ 4 , as seen in 
Fig. [10] 

In Trujillo Bueno & Fabiani Bendicho (1995), the 
Eddington approximation is used as the reference solution for a 
one -point angular quadrature n^-X with fi = ±1/ y3. The an- 
alytical expression of the Eddington approximation in a finite 
slab is: 



S E (e,r*,T)= l-(l-e) 



exp (- V3e t) + exp (- V3~6 (t* - r)J 
1 + Ve + (1 - Ve) exp (- V3eT") 
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Table 1. Comparison of results obtained with our (ALI+Ng) code and previous ones. Our N c is defined by e c = 0.01, while 
values from other authors are based on a graphical guess s c « 0.05. The optical thickness is t* = 2 x 10 8 . Note that in Trujillo 
Bueno & Manso Sainz (1999), N c values (in parenthesis) are given for a non-accelerated Jacobi scheme. These numbers have 
been divided by 2 in order to estimate the number of iterations when Ng acceleration is used. 
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(14) 

This is the exact solution of the monochromatic scattering 
problem when the mean intensity is calculated with the above- 
mentioned one-point angular quadrature. However, as is well- 
known, it gives only an approximation to the exact (i.e., multi- 
angle) solution of the full problem Q-Q- In other words, 
the true error given by Trujillo Bueno & Fabiani Bendicho 
(1995) is relative to the — 1 monochromatic scattering 
problem only, it does not give information on the error that 
would have been got by comparing the numerical solution to 
the tip = 1,3,5,... problem to the exact multi-angle solution 
(rtft = oo). In fact, as given in Table [2 for « A , = 1, when the 
solution of the n M — 1 problem is compared to the exact multi- 
angle solution, we find that the maximum error for n T — 9 is 
8.6 x 10~ 2 . The latter represents the maximum relative differ- 
ence between the Eddington approximation and the exact so- 
lution (Fig. II Q . A similar investigation, but for the two-level 
atom resonance-line scattering polarization problem, was car- 
ried out by Trujillo Bueno & Manso Sainz (1999), whose Table 
3 gives the surface true-error values of the fractional atomic po- 
larization for rip = 3,5,7, 11, . . ., 61. 




Fig. 11. Variation with t of the function #2(e,T*,T) = 
|1 - Sn(e, T*,r)/S(e, t*,t)| characterizing the accuracy of the 
Eddington approximation Se as given by J14i for the usual 
three couples (e, r*). 



6. Conclusion 

Our ALI code has been subjected to a wide range of tests, re- 
vealing at the same time its capabilities and its limits. Before 
developing these two points, we note that our conclusions are 
relative to the particular code we have used (based on Jacobi's 
method), specifically solving the standard problem Q-Q. The 
accuracy of the code is ultimately determined by the accuracy 
of the formal solver we have used (parabolic short characteris- 
tics). 

We have checked the great robustness of our code, which 
is certainly its most remarkable feature. It is able to solve the 
standard problem Q-Q for a wide range of input parameters 
e and r*, with no important lack of performance when e — > 
and/or t* — > +oo. 

However, the lowest accuracy of the ALI numerical solu- 
tions happens in the outermost layers of a star, corresponding 
to r lower than the thermalization depth 1 /k(e), these layers 
forming, by definition, the atmosphere of the star. The accu- 
racy of our code is not better than, say 10~ 2 , when we choose 
n T — 9, rtfj — 5 and limit the number of iterations to N < 100, 
as it is currently done in stellar atmospheres modelling. To im- 
prove the accuracy of the calculations up to ~ 10~ 3 , the param- 
eters n T , rtfj and N should be set to larger (but today impractical) 
values when solving the radiative transfer equation on a large 
frequency spectrum, i.e. at thousands of frequencies. Indeed 
we pointed out a truly noticeable improvement of the accu- 
racy when using finer grids in t or jj.. Such an observation was 
made easier by the use of a very accurate reference solution. Of 
course, increasing the level of refinement of both spatial and 
angular quadratures has a strong impact upon the number of 
iterations needed for convergence. However, to overcome this 
difficulty while keeping the same accuracy on the numerical 
solutions, methods based on Gauss-Seidel and successive over- 
relaxation iterations were already proposed by Trujillo Bueno 
& Fabiani Bendicho (1995). 

Another important question is relative to the propagation 
of errors in a stellar atmosphere model: to what extent are the 
main quantities provided by the model (populations of heavy 
particles, electron density, pressure, etc.) sensitive to the accu- 
racy on the RTE solution? We intend to tackle this subject in a 
future work by constructing an accurate - but still very ideal- 
ized - stellar atmosphere model, in which the main quantities 
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are first derived from an exact solution to the RTE, and then 
from the solution given by a ALI-based numerical method. 
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